//      ompr_main.cc
//      
//      Copyright 2012 tqlong <tqlong@espada>
//      
//      This program is free software; you can redistribute it and/or modify
//      it under the terms of the GNU General Public License as published by
//      the Free Software Foundation; either version 2 of the License, or
//      (at your option) any later version.
//      
//      This program is distributed in the hope that it will be useful,
//      but WITHOUT ANY WARRANTY; without even the implied warranty of
//      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//      GNU General Public License for more details.
//      
//      You should have received a copy of the GNU General Public License
//      along with this program; if not, write to the Free Software
//      Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
//      MA 02110-1301, USA.

#include <iostream>
#include <boost/program_options.hpp>
#include <boost/date_time/posix_time/posix_time.hpp> 

#include <armadillo>

#include "helper/stl_helper.h"
#include "helper/ublas_helper.h"
#include "helper/ml_data.h"

#include "ompr.h"

namespace po = boost::program_options;
po::options_description desc("Allowed options");
po::variables_map vm;

void process_options(int argc, char** argv);
void usage(int exitcode = 0);

void read_data();
void ompr();
void test();

#define NOW boost::posix_time::microsec_clock::local_time()

int main(int argc, char** argv)
{
  //MSG("Matching Pursuit for General Loss ... ");

  process_options(argc, argv);

  boost::posix_time::ptime time_start(boost::posix_time::microsec_clock::local_time());

  if (vm.count("help"))
    usage();
  else if (vm.count("read"))
    read_data();
  else if (vm.count("ompr"))
    ompr();
  else if (vm.count("test"))
    test();
  else
    usage();

  boost::posix_time::ptime time_end(boost::posix_time::microsec_clock::local_time());
  boost::posix_time::time_duration duration(time_end - time_start);
  MSG("Duration = " << duration);
}

void process_options(int argc, char** argv)
{
  desc.add_options()
    ("help", "produce this usage")
    ("input", po::value<std::string>(), "input file")
    ("output", po::value<std::string>(), "output file")
    ("debug", po::value<std::string>(), "debug file")
    ("seed", "random seed (REQUIRED)")
    ("fold", po::value<int>(), "number of folds")
    ("loss", po::value<std::string>(), "loss function (sq, log, smhinge, sqhinge")
    ("verbose", "print debug information")

    // OPTIMIZATION param
    ("optim_iter", po::value<int>(), "optim: maximum number of iterations")
    ("optim_tol", po::value<double>(), "optim: tolerance")
    ("optim_mem", po::value<int>(), "optim: L-BFGS memory")
    ("optim_eta", po::value<double>(), "optim: fixed step size")
    ("kernel", po::value<std::string>(), "Kernel cache maximum number of entries")
    ("kernel_cache", po::value<int>(), "Kernel cache maximum number of entries")

    // PARAMETERS
    ("mu", po::value<double>(), "the approximation parameter of smooth hinge loss")
    ("C", po::value<double>(), "loss penalty of SVM")
    ("last_label", "label at end of line")

    // OMPR parameter
    ("ompr_l", po::value<int>(), "OMPR number of replacement")
    ("ompr_k", po::value<int>(), "OMPR sparsity level")

    //JOB TYPE
    ("read", "read data file\n--input [--output]")
    ("ompr", "orthogonal matching pursuit with replacement\n"
             "--input [--output] [--ompr_l] [--ompr_k]")
    ("test", "test UBLAS")
    ;
  po::store(po::parse_command_line(argc, argv, desc), vm);
  po::notify(vm);
}

void usage(int exitcode)
{
  std::cout << desc << std::endl;
  exit(exitcode);
}


void read_data()
{
  
  std::string inputFile = vm["input"].as<std::string>();
  typedef DataTable<sparse_vec_t, sparse_matrix_t> Data;
  //typedef DataTable<vec_t, matrix_t> Data;

  Data D;
  DataReader::read_libsvm(inputFile, D, true); // add a bias term
  

  /*
  sparse_vec_t v(10);
  for (int i = 0; i < 100; i++) {
    if (i + 1 > sz(v)) v.resize(i+1);
    v.pb(i, i);
    MSG(sz(v) << " " << v.nnz());
  }
  */
  MSG("sz(D)=" << sz(D) << " dim(D)=" << D.dim());
  //MSG(D.data());
  //MSG(D.label());
}

void ompr() 
{
  std::string inputFile = vm["input"].as<std::string>();
  typedef DataTable<sparse_vec_t, sparse_matrix_t> Data;
  Data D;
  DataReader::read_libsvm(inputFile, D, false); // no bias term
  MSG(inputFile << ": sz(D)=" << sz(D) << " dim(D)=" << D.dim());
  OMPR<typeof(D)> solver(vm);
  solver.solve(D);
  
  const sparse_vec_t& x = solver.result();
  if (vm.count("output")) {
    std::string outFile = VM(vm, "output", std::string);
    FileAccessWrite(out, outFile);
    tr(x, it) out << it.index()+1 << ":" << *it << "\t "; // add 1 for libsvm convention
    out << ENDL; // out << x; it = x.begin() + 3
  }
  if (vm.count("verbose")) {
    tr(x, it) std::cerr << it.index()+1 << ":" << *it << "\t "; // add 1 for libsvm convention
    std::cerr << ENDL;
  }
  
}

struct node {
  int index;
  double val;
};
typedef std::vector<node> sp_t;

sp_t operator + (const sp_t& x, const sp_t& y) {
  sp_t v;
  
  int i = 0, j = 0;
  while (1) {
    if (j >= sz(y)) {
      v.insert(v.end(), x.begin()+i, x.end());
      break;
    }
    else if (i >= sz(x)) {
      v.insert(v.end(), y.begin()+j, y.end());
      break;
    }
    else if (x[i].index < y[j].index) {
      v.pb(x[i]);
      i++;
    }
    else if (x[i].index > y[j].index) {
      v.pb(y[j]);
      j++;
    }
    else {
      node a = x[i];
      a.val += y[j].val;
      v.pb(a);
      i++;
      j++;
    }
  }

  return v;
}

sp_t operator * (double s, const sp_t & x) {
  sp_t v = x;
  for (int i = 0; i < sz(x); i++) v[i].val = s*x[i].val;
  return v;
}

void test()
{
  int N = 100;

  boost::posix_time::ptime start = NOW;
  sparse_vec_t x(N), y(N);
  for (int i = 0; i < 2000000; i++) {
    y = x = scalar_vec_t(N, 1.1);
    x=x+7*y;
  }
  MSG("UBLAS Duration=" << NOW-start);

  start = NOW;
  sp_t xx(N), yy(N);
  for (int i = 0; i < 2000000; i++) {
    for (int j = 0; j < N; j++) xx[j].index = yy[j].index = j;
    for (int j = 0; j < N; j++) xx[j].val = yy[j].val = 1.1;

    xx = xx + 7 * yy;
  }
  MSG("HUA Duration=" << NOW-start);  

  start = NOW;
  vec_t xxx(N), yyy(N);
  for (int i = 0; i < 2000000; i++) {
    yyy = xxx = scalar_vec_t(N, 1.1);
    xxx=xxx+7*yyy;
  }
  MSG("UBLAS DENSE Duration=" << NOW-start);

  start = NOW;
  arma::sp_mat ax(N, 1), ay(N, 1);
  for (int i = 0; i < 2000000; i++) {
    for (int j = 0; j < N; j++) ay(j, 0) = ax(j, 0) = 1.1;
    ax = ax + 7 * ay;
  }
  MSG("ARMA SPARSE Duration=" << NOW-start);
  
}
